Validation#

In order to validate our results, this notebook compare our results with the flood mapping from the ‘Global Flood Monitoring’(GFM) for the same area and time range.

from importlib.resources import files

import holoviews as hv
import hvplot.xarray  # noqa
import numpy as np
import pystac_client
import xarray as xr
from dask.distributed import Client
from dask_flood_mapper import flood
from odc import stac as odc_stac
# define time range and area
time_range = "2023-10-11/2023-10-25"
bounding_box = [12.3, 54.3, 13.1, 54.6]
client = Client(processes=False, threads_per_worker=2, n_workers=1, memory_limit="12GB")

GFM#

# Connect to STAC catalog
eodc_catalog = pystac_client.Client.open("https://stac.eodc.eu/api/v1/")

# Search for available images

search = eodc_catalog.search(collections="GFM", bbox=bounding_box, datetime=time_range)
items_GFM = search.item_collection()

print(f"Found {len(items_GFM)} items")
Found 37 items

Load the data#

crs = "EPSG:4326"  # Coordinate Reference System - World Geodetic System 1984 (WGS84)
resolution = 0.00018  # 20 meter in degree

GFM_fd = odc_stac.load(
    items_GFM,
    bbox=bounding_box,
    crs=crs,
    bands=["tuw_likelihood", "tuw_flood_extent"],
    resolution=resolution,
    dtype="uint8",
    chunks={"x": 1000, "y": 1000, "time": -1},
)

Pre-processing#

# convert values 255 to NANs
GFM_fd["tuw_flood_extent"] = GFM_fd.tuw_flood_extent.where(
    GFM_fd.tuw_flood_extent != 255
).compute()

# select the band 'tuw_flood_extent'
GFM_fd = GFM_fd["tuw_flood_extent"]
GFM_fd
<xarray.DataArray 'tuw_flood_extent' (time: 28, latitude: 1668, longitude: 4445)> Size: 830MB
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
...
        ...,
        [nan, nan, nan, ...,  0.,  0.,  0.],
        [nan, nan, nan, ...,  0.,  0.,  0.],
        [nan, nan, nan, ...,  0.,  0.,  0.]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
Coordinates:
  * latitude     (latitude) float64 13kB 54.6 54.6 54.6 54.6 ... 54.3 54.3 54.3
  * longitude    (longitude) float64 36kB 12.3 12.3 12.3 12.3 ... 13.1 13.1 13.1
    spatial_ref  int32 4B 4326
  * time         (time) datetime64[ns] 224B 2023-10-11T05:33:18 ... 2023-10-2...
Attributes:
    nodata:   255

dask-flood-mapper#

# calculate flood decision
fd = flood.decision(bbox=bounding_box, datetime=time_range).compute()
fd
sigma naught datacube processed
harmonic parameter datacube processed
projected local incidence angle processed
<xarray.DataArray 'reproject-cd2f0e51db984733b1c696d947bda00b' (time: 9,
                                                                latitude: 1232,
                                                                longitude: 3283)> Size: 291MB
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [ 0.,  0.,  0., ..., nan, nan, nan],
        ...,
        [nan, nan,  0., ..., nan, nan, nan],
        [nan, nan,  0., ..., nan, nan, nan],
        [nan, nan,  0., ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ...,  0., nan, nan],
        ...,
        [nan, nan, nan, ...,  0.,  0.,  0.],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ...,  0., nan, nan],
        ...,
...
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [ 0.,  0.,  0., ..., nan, nan, nan],
        ...,
        [nan, nan,  0., ..., nan, nan, nan],
        [nan, nan,  0., ..., nan, nan, nan],
        [nan, nan,  0., ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ...,  0., nan, nan],
        ...,
        [nan, nan, nan, ...,  0.,  0.,  0.],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]])
Coordinates:
  * time         (time) datetime64[ns] 72B 2023-10-11T05:33:18 ... 2023-10-25...
  * latitude     (latitude) float64 10kB 54.6 54.6 54.6 54.6 ... 54.3 54.3 54.3
  * longitude    (longitude) float64 26kB 12.3 12.3 12.3 12.3 ... 13.1 13.1 13.1
    spatial_ref  int64 8B 0

Preprocessing#

Make both plots compatible to be plot together.

data_text = files("dask_flood_mapper.data").joinpath("wcover.tif")
wcover = xr.open_dataarray(data_text, band_as_variable=True).rename(
    {"y": "latitude", "x": "longitude"}
)
wcover
<xarray.DataArray 'band_1' (latitude: 1232, longitude: 3283)> Size: 16MB
[4044656 values with dtype=float32]
Coordinates:
  * longitude    (longitude) float64 26kB 12.3 12.3 12.3 12.3 ... 13.1 13.1 13.1
  * latitude     (latitude) float64 10kB 54.6 54.6 54.6 54.6 ... 54.3 54.3 54.3
    spatial_ref  int64 8B ...
Attributes:
    long_name:  wcover
fd = fd.where(wcover != 80)
fd
<xarray.DataArray 'reproject-cd2f0e51db984733b1c696d947bda00b' (time: 9,
                                                                latitude: 8,
                                                                longitude: 381)> Size: 219kB
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ...,  0.,  0.,  0.],
        [nan, nan, nan, ...,  0.,  0.,  0.],
        [nan, nan, nan, ...,  0.,  0.,  0.]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
...
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ...,  0.,  0.,  0.],
        [nan, nan, nan, ...,  0.,  0.,  0.],
        [nan, nan, nan, ...,  0.,  0.,  0.]]])
Coordinates:
  * time         (time) datetime64[ns] 72B 2023-10-11T05:33:18 ... 2023-10-25...
  * latitude     (latitude) float64 64B 54.34 54.33 54.32 ... 54.31 54.3 54.3
  * longitude    (longitude) float64 3kB 12.3 12.31 12.31 ... 13.1 13.1 13.1
    spatial_ref  int64 8B 0

Validation plot#

common_times = np.intersect1d(GFM_fd.time.values, fd.time.values)


def synced_plot(t):
    plot1 = GFM_fd.sel(time=t).hvplot.image(
        x="longitude",
        y="latitude",
        title="GFM flood map",
        cmap=["rgba(0, 0, 1, 0.1)", "darkblue"],
    )
    plot2 = fd.sel(time=t).hvplot.image(
        x="longitude",
        y="latitude",
        title="dask-flood-mapper",
        cmap=["rgba(0, 0, 1, 0.1)", "darkblue"],
    )

    return (plot1 + plot2).cols(2)


combined_plot = hv.DynamicMap(synced_plot, kdims="time").redim.values(time=common_times)

combined_plot
WARNING:param.Image00328: Image dimension(s) longitude and latitude are not evenly sampled to relative tolerance of 0.001. Please use the QuadMesh element for irregularly sampled data or set a higher tolerance on hv.config.image_rtol or the rtol parameter in the Image constructor.
WARNING:param.Image00328: Image dimension(s) longitude and latitude are not evenly sampled to relative tolerance of 0.001. Please use the QuadMesh element for irregularly sampled data or set a higher tolerance on hv.config.image_rtol or the rtol parameter in the Image constructor.
WARNING:param.Image00328: Image dimension(s) longitude and latitude are not evenly sampled to relative tolerance of 0.001. Please use the QuadMesh element for irregularly sampled data or set a higher tolerance on hv.config.image_rtol or the rtol parameter in the Image constructor.